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Abstract 

miRNAs are a class of small, single-stranded, non-coding RNAs that perform post-transcriptional repression of target 
genes by binding to 3' untranslated regions. Research has found that miRNAs involved in the regulation of many 
metabolic processes. Here we uncovered that the beef quality of Angus cattle sharply diversified after acute stress. By 
performing miRNA microarray analysis, 13 miRNAs were significantly differentially expressed in stressed group 
compared to control group. Using a bioinformatics method, 135 protein-coding genes were predicted as the targets of 
significant differentially expressed miRNAs. Gene Ontology (GO) term and Ingenuity Pathway Analysis (IPA) mined that 
these target genes involved in some important pathways, which may have impact on meat quality and beef tenderness. 
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Introduction 

MicroRNAs are one of the largest gene families and ac- 
count for -1% of the genome [1]. They are 21-25 nucleo- 
tide small, non-coding RNAs that post-transcriptionally 
repress the expression of protein-coding genes through 
binding to the 3' untranslated regions (UTR) of the target 
mRNAs [1-5]. Accumulated evidence indicates that miR- 
NAs are important in the regulation of many biological 
processes, such as developmental timing, cell metabolism, 
cell differentiation, cell death, cell proliferation, haemato- 
poiesis and patterning of the nervous system, etc [1,4,6]. 
Recent studies have uncovered muscle-specific miRNAs 
that regulate diverse aspects of muscle function, including 
myoblast proliferation, differentiation, contractility and 
stress responsiveness [7-10]. Disruption of miRNA biogen- 
esis causes diverse developmental defects, including ab- 
normal embryogenesis and depletion of stem cells [4]. It 
has been reported that micro RN A- 133a regulates cardio- 
myocyte proliferation and suppresses smooth muscle gene 
expression in heart [8]. miR-1 and miR-133 have distinct 
roles in modulating skeletal muscle proliferation and dif- 
ferentiation in cultured myoblast in vitro and in Xenopus 
laevis embryos in vivo [9]. miR-335 and miR-126 are 
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identified as metastasis suppressors in human breast can- 
cer because their expressions are lost in the majority of 
primary breast tumors [11]. Additionally, miRNAs have 
been found involved in viral infections, cancer, cardiovas- 
cular disease and neurological and muscular disorders 
[6,12-20]. With the progression of research, a large num- 
ber of miRNAs have been found to play roles in the 
regulation of metabolic process. Although there are 18226 
entries in miRBase, representing hairpin precursor miR- 
NAs and expressing 21643 mature miRNA products in 
168 species, only a handful of miRNAs have been studies 
deeply, and a range of functions extending beyond devel- 
opmental regulation have been revealed [4]. Especially, 
665 miRNAs in bovine are shown in the database, and 
some of them are studied in bovine cell in vitro, but few 
have been studied in vivo [21,22]. 

Beef tenderness is a complex characteristic influenced by 
many aspects, such as production, processing factors and 
cooking aspects, etc. More efforts have been focused on 
factors influencing meat quality, including breed, sex, feed, 
handling, environment, finishing weight and age at slaugh- 
ter, etc [23-28]. So far, no research is performed on whether 
the variation of beef tenderness is regulated by miRNAs. 
To test our hypothesis that acute stress may influence beef 
quality mediated by miRNAs, a miRNA microarray was 
used to detect differentially expressed miRNAs between 
stressed and non-stressed groups of cattle. The results from 
the study demonstrated that acute stress altered both beef 
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quality and miRNA expression, which will help us identify 
mechanisms underlying the control of beef tenderness. 

Results 

Differentially expressed miRNAs in LD muscle with 
differential stress status 

Warner-Bratzler shear forces (WBSF) measurements 
were made to evaluate variation of beef tenderness 
caused by acute stress. The results showed that the aver- 
age WBSF for stressed group and control group were 
19.74 kg and 5.04 kg, respectively. The stressed group 
was much tougher than the control (non-stressed) group 
from the student £-test result (P < 0.0001). To determine 
the miRNA expression patterns in Angus cattle with dif- 
ferent stress status, miRNA microarray analysis was con- 
ducted on LD muscle. These arrays were designed based 
on the miRBase Version 11.0 and contained 126 bovine 
miRNAs. For each miRNA, there were 4 to 8 repeat probes 
on each slide. After hybridization, washing, scanning, 
data were collected and then Limma package was ap- 
plied. A total of 13 miRNAs were identified as differen- 
tially expressed miRNAs using the criteria of P value less 
than 0.05 and FDR (false discover rate) less than 0.4 
(Table 1). Of these, one miRNA was down-regulated 
while 12 miRNAs were up-regulated in stressed group 
compared with control group. To reveal the overall ex- 
pression profiles of these differentially expressed miR- 
NAs in these two groups, clustering analysis was 
performed as previously described. The visualization 
showed that the expression patterns of these miRNAs 
can apparently divide these 6 individuals into stressed 
and control groups (Figure 1). To obtain a highly statisti- 
cally confident result, a stringent statistical significance 
threshold (P<0.05, Fold Change > 1.5) was used. Bta- 



Table 1 Differentially expressed miRNAs between 
stressed group and control in Angus cattle 



Name 


Fold Change 


P Value 


adj.P Value 


bta-miR-151 


0.839 


0.012 


0.397 


bta-miR-93 


1.128 


0.023 


0.397 


bta-miR-10b 


1.183 


0.037 


0.397 


bta-miR-20b 


1.187 


0.035 


0.397 


bta-let-7c 


1.200 


0.039 


0.397 


bta-let-7a 


1.206 


0.015 


0.397 


bta-miR-181b 


1.215 


0.043 


0.397 


bta-miR-99a 


1.226 


0.009 


0.397 


bta-miR-195 


1.244 


0.046 


0.397 


bta-miR-19b 


1.248 


0.027 


0.397 


bta-miR-660 


1.292 


0.020 


0.397 


bta-miR-125b 


1.321 


0.004 


0.397 


bta-miR-497 


1.624 


0.023 


0.397 
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Figure 1 Cluster analysis of significant differentially expressed 
miRNAs in microarray. These miRNAs were visualized with Treeview 
after hierarchical clustering. Each miRNA is represented by a single row 
of colored boxes; each individual from two groups is represented by a 
single column. Red color indicates up-regulated while green indicates 
down-regulated. miRNAs that were expressed at higher levels are 
assigned progressively brighter shades of red while miRNAs expressed 
at lower levels are assigned progressively brighter shades of green. 



Fold change was computed as (stressed/control). 



miR-497 was chosen as the most significantly differen- 
tially expressed miRNA to do further analysis. 

qPCR analysis of differentially expressed miRNA 

To validate the microarray results, mimic bta-miR-497 
were synthesized by miScript Primer Assays. Quantita- 
tive RT-PCR was performed to measure the expression 
level of bta-miR-497 in stressed and control groups. The 
results showed that the expression of bta-miR-497 
significantly increased in stressed group compared to 
control group (P < 0.05) (Figure 2), consistent with the 
miRNA microarray results (fold change = 1.62 in micro- 
array), namely, the expression of bta-miR-497 was 
increased after the acute stress. 

Prediction of targets of differentially expressed miRNA 
and function annotation 

To understand the potential functions of significantly 
differentially expressed miRNA in this diverse stress sta- 
tus, 135 genes were predicted as the potential targets of 
miR-497 in bovine by using bioinformatics method. To 
further explore the function of these predicted target 
genes, Gene Ontology analysis was performed. The 
results showed that the predicted target genes in GO 
biological process terms were enriched in cellular cata- 
bolic process and cellular process. In cellular component 
category, GO terms related to the cytoplasmic part, 
membrane-bounded organelle, intracellular membrane- 
bounded organelle, organelle, intracellular organelle, 
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Figure 2 Expression of bta-miR-497 in control and stressed groups 
was measured by qRT-PCR and normalized to U6. The quantitative 
results are represented as mean ±SEM (n = 3). A single asterisk {P <0.05). 



cytoplasm and intracellular organelle part. The molecular 
function category of GO terms showed that succinyltrans- 
ferase activity, purine nucleotide binding, ribonucleotide 
binding, purine ribonucleotide binding, purine ribonucleo- 
side triphosphate binding, GTP binding, guanyl nucleotide 
binding, guanyl ribonucleotide binding, S-acyltransf erase 
activity and GTPase activity were enriched. Summaries of 
the enriched GO term categories for predicted target 
genes are shown in the Table 2. 



To further visualize the pathways and networks these 
target genes related with, IPA of target genes was con- 
ducted. Analysis results showed that cell cycle, cell 
morphology, cellular function and maintenance, molecu- 
lar transport and cellular movement were ranked in 
the top of "Molecular and Cellular Functions". While, 
inhibition of angiogenesis by TSP1, D-glutamine and D- 
glutamate metabolism, G2/M DNA damage checkpoint 
regulation, galactose metabolism and nucleotide sugars 
metabolism were among the top canonical pathways. 
The most significant networks functioned in drug me- 
tabolism, endocrine system development and function, 
lipid metabolism, amino acid metabolism, molecular 
transport, small molecule biochemistry, gene expression, 
cellular movement, cell cycle, cardiovascular system de- 
velopment and function, organismal development, cancer 
and gastrointestinal disease. Summaries of the enriched 
networks and their functions are shown in Table 3 and 
graphical networks are represented (Figures 3, 4 and 5). 



Discussion 

Abnormal or disease conditions can induce dysregulation 
of mRNA and protein levels. It has been reported that 
muscle-specific miRNAs, miR-206 and miR-499, are 
upregulated and miR-1, miR-133a, and miR-133b are 
downregulated in extraocular muscles compared to limb 
muscle, concluding that a miRNA network contributes to 
the extraocular muscles by regulating posttranscriptional 



Table 2 Significant GO terms predicted target genes were involved in 



GOID 


Ontology 


Term 


q 


P 


GO:0044248 


biologicaLprocess 


cellular catabolic process 


13 


0.009486 


GO:0009987 


biologicaLprocess 


cellular process 


65 


0.043591 


GO:0044444 


cellular_component 


cytoplasmic part 


38 


0.008829 


GO:0043227 


cellular_component 


membrane-bounded organelle 


48 


0.03211 


GO:0043231 


cellular_component 


intracellular membrane- 
bounded organelle 


48 


0.03211 


GO:0043226 


cellular_component 


Organelle 


53 


0.035303 


GO:0043229 


cellular_component 


intracellular organelle 


53 


0.035303 


GO:0005737 


cellular_component 


Cytoplasm 


45 


0.041648 


GO:0044446 


cellular_component 


intracellular organelle part 


32 


0.044642 


GO:00 16748 


molecular_function 


succinyltransferase activity 


3 


0.00289 


GO:00 17076 


molecular_function 


purine nucleotide binding 


25 


0.009486 


GO:0032553 


molecular_function 


ribonucleotide binding 


25 


0.009486 


GO:0032555 


molecular_function 


purine ribonucleotide binding 


25 


0.009486 


GO:0035639 


molecular_function 


purine ribonucleoside 
triphosphate binding 


25 


0.009486 


GO:0005525 


molecular_function 


GTP binding 


11 


0.009486 


GO:00 19001 


molecular_function 


guanyl nucleotide binding 


11 


0.009876 


GO:0032561 


molecular_function 


guanyl ribonucleotide binding 


11 


0.009876 


GO:0016417 


molecular_function 


S-acyltransferase activity 


3 


0.013341 


GO:0003924 


molecular_function 


GTPase activity 


7 


0.03211 
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Table 3 Networks and functions that target genes are related with 



ID 


Score 


Focus Molecules 


Top Functions 


1 


19 




Drug Metabolism, Endocrine System Development and 
Function, Lipid Metabolism 


2 


15 


11 


Amino Acid Metabolism, Molecular Transport, Small 
Molecule Biochemistry 


3 


15 


11 


Gene Expression, Cellular Movement, Cell Cycle 


4 


15 


11 


Cardiovascular System Development and Function, 
Organismal Development, Cancer 


5 


15 


11 


Cell Cycle, Cancer, Gastrointestinal Disease 


6 


11 




Cardiovascular Disease, Gene Expression, 
Hematological Disease 


7 


2 


1 


Carbohydrate Metabolism, Lipid Metabolism, 
Small Molecule Biochemistry 


8 


2 


i 


Cellular Growth and Proliferation, Gene Expression, 
Infectious Disease 


9 


2 


i 


RNA Post-Transcriptional Modification, Cellular 
Compromise, Cellular Development 


10 


2 


1 


Cancer, Embryonic Development, Neurological Disease 


11 


2 


1 


RNA Damage and Repair, Nutritional Disease, 
Organismal Injury and Abnormalities 


12 


2 


1 


Cellular Growth and Proliferation, Developmental Disorder, 
Embryonic Development 


13 


2 


1 


Developmental Disorder, Genetic Disorder, Metabolic Disease 


14 


2 


1 


Cell Signaling, Cellular Assembly and Organization, 
Cellular Function and Maintenance 


15 


2 




Post-Translational Modification, Protein Synthesis, Cell-To-Cell 
Signaling and Interaction 


16 


2 




Cellular Assembly and Organization, Cellular Function and 
Maintenance, Cellular Movement 


17 


2 




Cellular Assembly and Organization, Cellular Function and 
Maintenance, Cell-To-Cell Signaling and Interaction 


18 


2 




Cell Cycle, Cellular Movement, Embryonic Development 


19 


2 




Infectious Disease, DNA Replication, Recombination, and Repair, 
Gene Expression 



ID represents the network No. Score means gene number in this network and Focus molecular means the target gene number in this network. 



expression of genes involved in structure, signaling, me- 
tabolism, angiogenesis, myogenesis, and regeneration in 
extraocular muscles [7]. In addition, miR-145 is found to 
be necessary for myocardin-induced reprogramming of 
adult fibroblasts into smooth muscle cells and can induce 
differentiation of multipotent neural crest stem cells into 
vascular smooth muscle [10]. Meanwhile, miR-145 and 
miR-143 cooperatively target a network of transcription 
factors to promote differentiation and repress proliferation 
of smooth muscle cells [10]. Both also act as integral com- 
ponents of the regulatory network whereby serum re- 
sponse factor controls cytoskeletal remodeling and 
phenotypic switching of smooth muscle cells during vas- 
cular disease [29]. In our study, several miRNAs were 
found to be dysregulated due to different stress status, of 
which, some have been previously studied. For example, 
miR-497 has been found to promote ischemic neuronal 
death by negatively regulating antiapoptotic proteins [30]. 



Another research found that miR-497 and miR-302b co- 
regulate ethanol-induced neuronal cell death through 
BCL2 protein and cyclin D2 [31]. But its function in 
muscle development has not been reported yet. Therefore, 
these finding further suggest that miRNAs may play some 
roles on transcriptional circuits controlling gene expres- 
sion in skeletal muscle. 

Notably, the surgical implantation of rumen canulas 
imitated a non fatal form of hardware disease. Hardware 
disease occurs when an animal ingests a sharp piece of 
metal and the metal pierces the rumen or reticulum wall. 
As expected, the phenotype in this study indicated that 
those animals undergoing this stress had significantly 
higher WBSF. In this research, we identified differently 
expressed miRNAs associated with divergent stress sta- 
tus in LD muscles samples between stressed and control 
groups. The annotation of predicted target genes further 
showed that miRNA may be involved in important 
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pathways regulating target genes, such as lipid metabo- 
lism, amino acid metabolism, gene expression, molecular 
transport, etc. In the future, the predicted miRNA targets 
need to be validated in vitro and the expression levels of 
corresponding target genes and proteins should be mea- 
sured, which will help to elucidate how miRNAs regulate 
gene transcription and protein expression in the vari- 
ation of beef quality and tenderness. 

Materials and methods 

Sample collection and experiment design 

Seven purebred Angus steers were obtained from Wye 
Angus farm (Queenstown, MD). After weaning the steers 
were acclimated to a pelleted forage diet only to meet 
maintenance needs. At 10 months of age, 4 steers under- 
went a surgical procedure that involved anesthetization 
and placement of a rumen catheter. The surgery was 
acute stress compared to normal growth condition. 
Three steers that received no surgery were designated as 
control group. At the age of 1 year, the steers were har- 
vested. After harvest 10 mg longissimus dorsi (LD) 
muscle from the 12 th to 13 th rib of the right side of the 
carcass were placed in RNAlater solution (Qiagen, Valen- 
cia, CA) and stored at -80°C for further analysis. Steaks 



of the LD from the 12 th to 13 th rib of the left side of the 
carcass were obtained, vacuum packed, stored at 4°C for 
a total of 14 days post harvest, and then frozen at -20°C. 
Once all steaks were obtained, aged, the steaks were 
thawed at 4°C, cooked to an internal temperature at 70°C, 
cooled, cored and then analyzed for WBSF as previously 
described [32]. After WBSF data were analyzed by 
student £-test, three extremely tough individuals were 
chosen to be designed as stressed group and three cattle 
without stress were designed as control group. Based on 
these tough and control groups, a total of 6 miRNA 
microarrays were hybridized and analyzed. All proce- 
dures were approved by the University of Maryland Insti- 
tutional Animal Care and Use Committee (Protocol # R- 
07-05). 

RNA extraction and miRNA microarray hybridization 

Total RNAs from the 6 samples were extracted using 
miRNAeasy Mini Kit (Qiagen) as described in the manu- 
facturer s instructions. The RNAs were quantified by Nano- 
Drop ND 1000 Spectrophotometer (Thermo-scientific, 
Wilmington, DE) and RNA integrity determined by 2100 
Bioanalyzer (Agilent Technologies, Santa Clara, CA). Equal 
aliquots of total RNA from each sample were pooled 
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Figure 4 The top 2# network target genes involved. Solid line represents direct interaction and dash line represents indirect interaction. 



together as common reference RNA. One (ig total RNA 
from each sample or common reference were labeled with 
Hy3™ and Hy5™ fluorescent label, respectively, with the 
help of the miRCURY™ LAN Array power labeling kit 
(Exiqon, Denmark) following the instructions. The Hy3™- 
labeled samples and a Hy5™-labeled reference RNA 
sample were mixed pair-wise and hybridized to the miR- 
CURY™ LNA array (Version 9.2; Exiqon, Denmark), 
which contained capture probes targeting all of the miR- 
NAs for all the species registered in the miRBase (Version 
11.0) at the Sanger Institute. One hundred and twenty-six 
of these probes are bovine-related miRNAs in the miRBase 
version. Hybridization was performed according to the 
miRCURY™ LNA array manual with a Tecan HS4800 
hybridization station (Tecan, Austria). After hybridization, 
the microarray slides were scanned and stored in an ozone 
free environment to prevent potential bleaching of the 
fluorescent dyes. The miRCURY™ LNA array microarray 
slides were scanned using the Agilent G2565BA Microarray 
Scanner System (Agilent) and image analysis was performed 
with ImaGene 8.0 software (BioDiscovery, Inc., USA). 

miRNA microarray data analysis 

Microarray data were analyzed in R using the Linear 
models for microarray data (Limma) package. For each 
miRNA, quantified signals within arrays were averaged. 
Normalizations within arrays and between arrays were 



performed using the global LOWESS (LOcally WEighted 
Scatterplot Smoothing) regression algorithm. Contrasts 
were made to compare stressed and control groups. Dif- 
ferentially expressed miRNAs were selected to do further 
analysis using the stringent statistic criteria of p value 
less than 0.05 and FDR (false discover rate) less than 0.4. 

qRT-PCR analysis of miRNA expression 

Total mRNAs including miRNAs were extracted from 6 
same samples using miRNeasy Mini Kit (QIAGEN) and 
RNeasy Mini Kit (QIAGEN) according to the standard 
protocol. mRNAs were reversely transcribed and quanti- 
fied with miScript Reverse Transcription Kit (QIAGEN), 
miScript SYBR Green PCR Kit (QIAGEN), and miScript 
Primer assays (QIAGEN). In the reverse transcription 
control, PCR water (Invitrogen) was used to replace 
miRNA samples. Briefly, l(ig of purified miRNA was 
used for reverse transcription, and then diluted to 5 
volumes. Two \i\ of diluted RT products were used for 
real-time PCR quantification. Two types of controls were 
applied in real-time PCR, including reverse transcription 
control and blank using PCR water, to ensure that no 
amplicon was observed in the controls. U6 were used 
as normalization controls. Data were analyzed using the 
2-AACT met ] loc [ an d student T tests were used to com- 
pare the miRNAs expression levels (SAS version 9.2). 
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Figure 5 The top 3# network target genes involved. Solid line represents direct interaction and dash line represents indirect interaction. 



Here we only validated the most significant miRNA, 
namely bta-miR-497, which sequence is shown as CAG- 
CAGCACACUGUGGUUUGUA. The mimic miRNA of 
bta-miR-497 was synthesized by Qiagen. 

Prediction of miRNA targets 

The target genes for miRNAs were predicted by Target- 
ScanHuman (http://www.targetscan.org/vert_50/). In the 
menu of "Select a species", cow was chosen and the 
names of significantly differentially expressed miRNAs 
were inputted and then submitted. From the output only 
the genes with the conserved sites were reserved as pre- 
dicted target genes of this miRNA. 

Data mining and network analysis of significantly 
differentially expressed miRNAs and predicted target 
genes 

Hierarchical clustering of significantly differentially ex- 
pressed miRNAs was performed using Cluster 3.0 [33]. 
The expression data were further filtered, adjusted and 
normalized. Average linkage clustering was performed 
and visualized using Treeview. The initial information on 
Gene Ontology [15] functions and functional relevance 
of predicted target genes was obtained from Gene Ontol- 
ogy Enrichment Analysis Software Toolkit (GOEAST) 



[34]. The GO analysis included biological process, mo- 
lecular function and cellular component. Ingenuity Path- 
way Analysis (IPA, Ingenuity System, Redwood City, CA) 
was used to generate networks and assess statistically 
relevant biofunctions and canonical pathways that pre- 
dicted target genes are involved in. These genes were 
mapped to corresponding genes in the Ingenuity know- 
ledge database. The biofunctional analysis identified the 
molecular and cellular function, physiological system de- 
velopment and function. Canonical Pathway Analysis 
identified the most significant pathways in the dataset. 
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